Discriminant Analysis

Linear Discriminant Analysis (LDA) y Quadratic Discriminant Analysis (QDA)

Author

Dante Conti, Sergi Ramirez, (c) IDEAI

1 Definición del problema

1.1 Contexto

En este documento trabajaremos con una base de datos clásica de clasificación: Wine Recognition Dataset. Cada fila representa una muestra de vino procedente de una de tres variedades distintas. El objetivo es construir modelos capaces de clasificar correctamente la variedad del vino a partir de medidas químicas.

El fichero original LDA.R desarrolla un flujo completo de Linear Discriminant Analysis (LDA) en R. En este documento se mantiene esa lógica, se amplía la explicación docente, se organiza el material siguiendo la estructura del ejemplo KnnNB.qmd y se añade una sección completa de Quadratic Discriminant Analysis (QDA).

Las variables de la base de datos son:

Wine: clase o variedad del vino. Es la variable respuesta.
Alcohol: contenido de alcohol.
Malic.acid: ácido málico.
Ash: cenizas.
Acl: alcalinidad de las cenizas.
Mg: magnesio.
Phenols: fenoles totales.
Flavanoids: flavonoides.
Nonflavanoid.phenols: fenoles no flavonoides.
Proanth: proantocianinas.
Color.int: intensidad del color.
Hue: tonalidad.
OD: medida de dilución óptica.
Proline: prolina.

1.2 Objetivo

El objetivo es enseñar cómo funcionan LDA y QDA desde un punto de vista aplicado:

Cargar y preparar una base de datos de clasificación multiclase.
Dividir los datos en entrenamiento y test.
Ajustar un modelo LDA.
Interpretar las funciones discriminantes lineales.
Visualizar las clases en el espacio discriminante.
Clasificar observaciones nuevas.
Evaluar el modelo con matriz de confusión y accuracy.
Ajustar un modelo QDA.
Comparar LDA y QDA desde el punto de vista conceptual y predictivo.

2 Fundamentos teóricos

2.1 Idea general del análisis discriminante

El análisis discriminante es una familia de métodos supervisados de clasificación. Partimos de observaciones etiquetadas, es decir, observaciones para las que conocemos la clase real. A partir de esas observaciones, el modelo aprende una regla de clasificación que después se usa para asignar clase a nuevas observaciones.

En este caso queremos predecir:

\[ Y \in \{1, 2, 3\} \]

usando un vector de variables predictoras:

\[ X = (X_1, X_2, \ldots, X_p) \]

La pregunta central es:

Dadas las características químicas de un vino, ¿a qué clase pertenece?

2.2 LDA: Linear Discriminant Analysis

LDA asume que las variables predictoras siguen una distribución normal multivariante dentro de cada grupo y que todos los grupos comparten la misma matriz de covarianzas.

De forma simplificada:

\[ X \mid Y = k \sim \mathcal{N}(\mu_k, \Sigma) \]

Donde:

$\mu_k$ es el vector de medias del grupo $k$.
$\Sigma$ es la matriz de covarianzas común para todos los grupos.

La consecuencia práctica es que las fronteras de decisión entre clases son lineales. Por eso el método se llama análisis discriminante lineal.

2.3 QDA: Quadratic Discriminant Analysis

QDA relaja una de las hipótesis fuertes de LDA. En QDA se permite que cada clase tenga su propia matriz de covarianzas:

\[ X \mid Y = k \sim \mathcal{N}(\mu_k, \Sigma_k) \]

Donde:

$\mu_k$ es el vector de medias del grupo $k$.
$\Sigma_k$ es la matriz de covarianzas específica del grupo $k$.

La consecuencia práctica es que las fronteras de decisión pueden ser cuadráticas. QDA es más flexible que LDA, pero también necesita más datos para estimar correctamente más parámetros.

2.4 Diferencias entre LDA y QDA

Aspecto	LDA	QDA
Covarianzas	Comunes para todas las clases	Una covarianza por clase
Frontera de decisión	Lineal	Cuadrática
Flexibilidad	Menor	Mayor
Riesgo de sobreajuste	Menor	Mayor
Requisitos de datos	Menos exigente	Más exigente
Interpretabilidad	Alta	Media

En general:

Si las clases tienen estructuras de dispersión parecidas, LDA suele ser una buena opción.
Si las clases tienen estructuras de dispersión muy distintas y hay suficientes datos, QDA puede capturar mejor la separación.

3 Carga y preparación de datos

La base de datos se carga desde una URL pública. Para que el documento sea reproducible, todo el proceso de carga, transformación y partición se incluye dentro del propio .qmd.

R
Python

Mostrar código

# ==============================================================================
# Librerías necesarias
# ==============================================================================

# MASS contiene las funciones lda() y qda().
library(MASS)

# caret se utilizará para crear la partición train/test y calcular matrices
# de confusión de forma cómoda.
library(caret)

# ggplot2 se usará para visualizaciones más docentes.
library(ggplot2)

# dplyr se usará para manipulación sencilla de datos.
library(dplyr)

# klaR permite dibujar regiones de decisión con partimat().
# Si no está instalado, se omite la parte correspondiente más adelante.
if (requireNamespace("klaR", quietly = TRUE)) {
  library(klaR)
}

# biotools permite realizar el test de Box's M para comparar matrices de
# covarianzas entre grupos. Si no está instalado, se omite esa parte.
if (requireNamespace("biotools", quietly = TRUE)) {
  library(biotools)
}

Mostrar código

# ==============================================================================
# Carga de la base de datos Wine
# ==============================================================================

url <- "https://gist.githubusercontent.com/tijptjik/9408623/raw/b237fa5848349a14a14e5d4107dc7897c21951f5/wine.csv"

wine <- read.csv(url, sep = ",", stringsAsFactors = TRUE)

# Comprobamos los nombres de las variables.
names(wine)

 [1] "Wine"                 "Alcohol"              "Malic.acid"          
 [4] "Ash"                  "Acl"                  "Mg"                  
 [7] "Phenols"              "Flavanoids"           "Nonflavanoid.phenols"
[10] "Proanth"              "Color.int"            "Hue"                 
[13] "OD"                   "Proline"

Mostrar código

# La variable Wine es la clase del problema. Nos aseguramos de que sea factor.
wine$Wine <- as.factor(wine$Wine)

# Vista inicial de la base de datos.
str(wine)

'data.frame':   178 obs. of  14 variables:
 $ Wine                : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
 $ Alcohol             : num  14.2 13.2 13.2 14.4 13.2 ...
 $ Malic.acid          : num  1.71 1.78 2.36 1.95 2.59 1.76 1.87 2.15 1.64 1.35 ...
 $ Ash                 : num  2.43 2.14 2.67 2.5 2.87 2.45 2.45 2.61 2.17 2.27 ...
 $ Acl                 : num  15.6 11.2 18.6 16.8 21 15.2 14.6 17.6 14 16 ...
 $ Mg                  : int  127 100 101 113 118 112 96 121 97 98 ...
 $ Phenols             : num  2.8 2.65 2.8 3.85 2.8 3.27 2.5 2.6 2.8 2.98 ...
 $ Flavanoids          : num  3.06 2.76 3.24 3.49 2.69 3.39 2.52 2.51 2.98 3.15 ...
 $ Nonflavanoid.phenols: num  0.28 0.26 0.3 0.24 0.39 0.34 0.3 0.31 0.29 0.22 ...
 $ Proanth             : num  2.29 1.28 2.81 2.18 1.82 1.97 1.98 1.25 1.98 1.85 ...
 $ Color.int           : num  5.64 4.38 5.68 7.8 4.32 6.75 5.25 5.05 5.2 7.22 ...
 $ Hue                 : num  1.04 1.05 1.03 0.86 1.04 1.05 1.02 1.06 1.08 1.01 ...
 $ OD                  : num  3.92 3.4 3.17 3.45 2.93 2.85 3.58 3.58 2.85 3.55 ...
 $ Proline             : int  1065 1050 1185 1480 735 1450 1290 1295 1045 1045 ...

Mostrar código

summary(wine)

 Wine      Alcohol        Malic.acid         Ash             Acl       
 1:59   Min.   :11.03   Min.   :0.740   Min.   :1.360   Min.   :10.60  
 2:71   1st Qu.:12.36   1st Qu.:1.603   1st Qu.:2.210   1st Qu.:17.20  
 3:48   Median :13.05   Median :1.865   Median :2.360   Median :19.50  
        Mean   :13.00   Mean   :2.336   Mean   :2.367   Mean   :19.49  
        3rd Qu.:13.68   3rd Qu.:3.083   3rd Qu.:2.558   3rd Qu.:21.50  
        Max.   :14.83   Max.   :5.800   Max.   :3.230   Max.   :30.00  
       Mg            Phenols        Flavanoids    Nonflavanoid.phenols
 Min.   : 70.00   Min.   :0.980   Min.   :0.340   Min.   :0.1300      
 1st Qu.: 88.00   1st Qu.:1.742   1st Qu.:1.205   1st Qu.:0.2700      
 Median : 98.00   Median :2.355   Median :2.135   Median :0.3400      
 Mean   : 99.74   Mean   :2.295   Mean   :2.029   Mean   :0.3619      
 3rd Qu.:107.00   3rd Qu.:2.800   3rd Qu.:2.875   3rd Qu.:0.4375      
 Max.   :162.00   Max.   :3.880   Max.   :5.080   Max.   :0.6600      
    Proanth        Color.int           Hue               OD       
 Min.   :0.410   Min.   : 1.280   Min.   :0.4800   Min.   :1.270  
 1st Qu.:1.250   1st Qu.: 3.220   1st Qu.:0.7825   1st Qu.:1.938  
 Median :1.555   Median : 4.690   Median :0.9650   Median :2.780  
 Mean   :1.591   Mean   : 5.058   Mean   :0.9574   Mean   :2.612  
 3rd Qu.:1.950   3rd Qu.: 6.200   3rd Qu.:1.1200   3rd Qu.:3.170  
 Max.   :3.580   Max.   :13.000   Max.   :1.7100   Max.   :4.000  
    Proline      
 Min.   : 278.0  
 1st Qu.: 500.5  
 Median : 673.5  
 Mean   : 746.9  
 3rd Qu.: 985.0  
 Max.   :1680.0

Mostrar código

head(wine)

Mostrar código

# ==============================================================================
# Carga de la base de datos Wine en Python
# ==============================================================================

import pandas as pd
import numpy as np

url = "https://gist.githubusercontent.com/tijptjik/9408623/raw/b237fa5848349a14a14e5d4107dc7897c21951f5/wine.csv"

wine_py = pd.read_csv(url)

# Convertimos la variable objetivo a tipo categórico.
wine_py["Wine"] = wine_py["Wine"].astype("category")

print(wine_py.head())

  Wine  Alcohol  Malic.acid   Ash  ...  Color.int   Hue    OD  Proline
0    1    14.23        1.71  2.43  ...       5.64  1.04  3.92     1065
1    1    13.20        1.78  2.14  ...       4.38  1.05  3.40     1050
2    1    13.16        2.36  2.67  ...       5.68  1.03  3.17     1185
3    1    14.37        1.95  2.50  ...       7.80  0.86  3.45     1480
4    1    13.24        2.59  2.87  ...       4.32  1.04  2.93      735

[5 rows x 14 columns]

Mostrar código

print(wine_py.info())

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 178 entries, 0 to 177
Data columns (total 14 columns):
 #   Column                Non-Null Count  Dtype   
---  ------                --------------  -----   
 0   Wine                  178 non-null    category
 1   Alcohol               178 non-null    float64 
 2   Malic.acid            178 non-null    float64 
 3   Ash                   178 non-null    float64 
 4   Acl                   178 non-null    float64 
 5   Mg                    178 non-null    int64   
 6   Phenols               178 non-null    float64 
 7   Flavanoids            178 non-null    float64 
 8   Nonflavanoid.phenols  178 non-null    float64 
 9   Proanth               178 non-null    float64 
 10  Color.int             178 non-null    float64 
 11  Hue                   178 non-null    float64 
 12  OD                    178 non-null    float64 
 13  Proline               178 non-null    int64   
dtypes: category(1), float64(11), int64(2)
memory usage: 18.5 KB
None

4 Análisis exploratorio inicial

Antes de ajustar cualquier modelo conviene inspeccionar si las clases están razonablemente separadas en algunas variables y si hay diferencias claras en las medias.

R
Python

Mostrar código

# Número de observaciones por clase.
table(wine$Wine)


 1  2  3 
59 71 48

Mostrar código

# Medias por grupo para algunas variables relevantes.
wine %>%
  group_by(Wine) %>%
  summarise(
    Alcohol_media = mean(Alcohol),
    Flavanoids_media = mean(Flavanoids),
    Color_int_media = mean(Color.int),
    Proline_media = mean(Proline),
    .groups = "drop"
  )

Mostrar código

# Boxplot de Alcohol por clase.
ggplot(wine, aes(x = Wine, y = Alcohol, fill = Wine)) +
  geom_boxplot(alpha = 0.75) +
  theme_minimal() +
  labs(
    title = "Distribución de Alcohol por clase de vino",
    x = "Clase de vino",
    y = "Alcohol"
  ) +
  theme(legend.position = "none")

Mostrar código

# Relación entre flavonoides y prolina.
ggplot(wine, aes(x = Flavanoids, y = Proline, color = Wine)) +
  geom_point(size = 2, alpha = 0.85) +
  theme_minimal() +
  labs(
    title = "Separación visual de clases usando dos variables",
    x = "Flavanoids",
    y = "Proline"
  )

Mostrar código

# Resumen de frecuencias por clase.
print(wine_py["Wine"].value_counts())

Wine
2    71
1    59
3    48
Name: count, dtype: int64

Mostrar código

# Medias por grupo.
print(
    wine_py.groupby("Wine", observed=True)[["Alcohol", "Flavanoids", "Color.int", "Proline"]]
    .mean()
)

        Alcohol  Flavanoids  Color.int      Proline
Wine                                               
1     13.744746    2.982373   5.528305  1115.711864
2     12.278732    2.080845   3.086620   519.507042
3     13.153750    0.781458   7.396250   629.895833

5 Partición en entrenamiento y test

La partición train/test sirve para estimar el rendimiento del modelo sobre observaciones que no se han utilizado durante el ajuste. En este documento usaremos una partición estratificada para conservar la proporción de clases.

R
Python

Mostrar código

# ==============================================================================
# Partición train/test
# ==============================================================================

set.seed(123)

# createDataPartition mantiene aproximadamente la proporción de cada clase.
trainIndex <- createDataPartition(wine$Wine, p = 0.8, list = FALSE, times = 1)

train <- wine[trainIndex, ]
test  <- wine[-trainIndex, ]

# Comprobamos el tamaño de cada conjunto.
dim(train)

[1] 144  14

Mostrar código

dim(test)

[1] 34 14

Mostrar código

# Comprobamos la distribución de clases.
prop.table(table(train$Wine))


        1         2         3 
0.3333333 0.3958333 0.2708333

Mostrar código

prop.table(table(test$Wine))


        1         2         3 
0.3235294 0.4117647 0.2647059

Mostrar código

from sklearn.model_selection import train_test_split

X_py = wine_py.drop(columns=["Wine"])
y_py = wine_py["Wine"]

X_train_py, X_test_py, y_train_py, y_test_py = train_test_split(
    X_py,
    y_py,
    test_size=0.2,
    random_state=123,
    stratify=y_py
)

print(X_train_py.shape)

(142, 13)

Mostrar código

print(X_test_py.shape)

(36, 13)

Mostrar código

print(y_train_py.value_counts(normalize=True))

Wine
2    0.401408
1    0.330986
3    0.267606
Name: proportion, dtype: float64

Mostrar código

print(y_test_py.value_counts(normalize=True))

Wine
2    0.388889
1    0.333333
3    0.277778
Name: proportion, dtype: float64

6 LDA

6.1 Ajuste del modelo LDA

El modelo LDA se ajusta con la función lda() del paquete MASS. La fórmula Wine ~ . indica que queremos predecir Wine usando todas las demás variables como predictoras.

R
Python

Mostrar código

# ==============================================================================
# Ajuste del modelo LDA
# ==============================================================================

wine_lda <- lda(Wine ~ ., data = train)

# Mostramos el objeto del modelo.
wine_lda

Call:
lda(Wine ~ ., data = train)

Prior probabilities of groups:
        1         2         3 
0.3333333 0.3958333 0.2708333 

Group means:
   Alcohol Malic.acid      Ash      Acl        Mg  Phenols Flavanoids
1 13.75771   2.022917 2.457083 16.82708 105.39583 2.833542  2.9783333
2 12.28614   1.914211 2.247368 20.44035  92.35088 2.272105  2.1182456
3 13.17974   3.434615 2.450769 21.44872  98.25641 1.704615  0.7966667
  Nonflavanoid.phenols  Proanth Color.int       Hue       OD   Proline
1            0.2872917 1.904583  5.455000 1.0562500 3.171458 1115.3542
2            0.3721053 1.604035  3.138070 1.0497544 2.775439  512.5965
3            0.4441026 1.149487  7.352308 0.6815385 1.701795  618.8462

Coefficients of linear discriminants:
                              LD1           LD2
Alcohol              -0.471728928  0.7636107428
Malic.acid            0.237567323  0.4210833071
Ash                  -0.522618904  2.5552627175
Acl                   0.178082167 -0.1514517232
Mg                   -0.016765890  0.0002713628
Phenols               0.747281803  0.0819012804
Flavanoids           -1.610806838 -0.6869403889
Nonflavanoid.phenols -1.615495361 -1.7760864578
Proanth              -0.176937781 -0.1745316895
Color.int             0.436276775  0.2745286415
Hue                  -0.120521287 -1.0831143233
OD                   -0.951293724  0.0856652866
Proline              -0.003160704  0.0026216255

Proportion of trace:
   LD1    LD2 
0.6971 0.3029

Mostrar código

# ==============================================================================
# Elementos principales del objeto LDA
# ==============================================================================

# 1) Probabilidades a priori: proporción de cada clase en train.
wine_lda$prior

        1         2         3 
0.3333333 0.3958333 0.2708333

Mostrar código

# 2) Medias de grupo: centro de cada clase para cada variable.
wine_lda$means

   Alcohol Malic.acid      Ash      Acl        Mg  Phenols Flavanoids
1 13.75771   2.022917 2.457083 16.82708 105.39583 2.833542  2.9783333
2 12.28614   1.914211 2.247368 20.44035  92.35088 2.272105  2.1182456
3 13.17974   3.434615 2.450769 21.44872  98.25641 1.704615  0.7966667
  Nonflavanoid.phenols  Proanth Color.int       Hue       OD   Proline
1            0.2872917 1.904583  5.455000 1.0562500 3.171458 1115.3542
2            0.3721053 1.604035  3.138070 1.0497544 2.775439  512.5965
3            0.4441026 1.149487  7.352308 0.6815385 1.701795  618.8462

Mostrar código

# 3) Coeficientes de las funciones discriminantes.
# En un problema con K clases, como máximo hay K - 1 funciones discriminantes.
wine_lda$scaling

                              LD1           LD2
Alcohol              -0.471728928  0.7636107428
Malic.acid            0.237567323  0.4210833071
Ash                  -0.522618904  2.5552627175
Acl                   0.178082167 -0.1514517232
Mg                   -0.016765890  0.0002713628
Phenols               0.747281803  0.0819012804
Flavanoids           -1.610806838 -0.6869403889
Nonflavanoid.phenols -1.615495361 -1.7760864578
Proanth              -0.176937781 -0.1745316895
Color.int             0.436276775  0.2745286415
Hue                  -0.120521287 -1.0831143233
OD                   -0.951293724  0.0856652866
Proline              -0.003160704  0.0026216255

Mostrar código

# 4) Proporción de traza: importancia relativa de cada función discriminante.
wine_lda$svd^2 / sum(wine_lda$svd^2)

[1] 0.6971412 0.3028588

Mostrar código

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

# ==============================================================================
# Ajuste del modelo LDA en Python
# ==============================================================================

lda_py = LinearDiscriminantAnalysis()
lda_py.fit(X_train_py, y_train_py)

LinearDiscriminantAnalysis()

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Mostrar código

print("Clases:", lda_py.classes_)

Clases: [1 2 3]

Mostrar código

print("Probabilidades a priori:", lda_py.priors_)

Probabilidades a priori: [0.33098592 0.40140845 0.26760563]

Mostrar código

print("Proporción explicada por cada eje discriminante:", lda_py.explained_variance_ratio_)

Proporción explicada por cada eje discriminante: [0.63661093 0.36338907]

6.2 Proyección en el espacio discriminante

LDA no solo clasifica. También construye nuevas variables llamadas funciones discriminantes. Estas funciones son combinaciones lineales de las variables originales y buscan maximizar la separación entre clases.

Como hay tres clases, LDA puede construir como máximo dos funciones discriminantes:

\[ K - 1 = 3 - 1 = 2 \]

R
Python

Mostrar código

# ==============================================================================
# Proyección de train sobre los ejes discriminantes
# ==============================================================================

lda_train_pred <- predict(wine_lda, train[, -1])

# lda_train_pred$x contiene las coordenadas LD1 y LD2.
head(lda_train_pred$x)

        LD1       LD2
1 -5.200632 1.7904175
2 -4.304108 0.9163969
3 -3.758153 1.4211415
5 -1.731594 0.5375356
6 -4.941377 2.8688266
7 -4.794154 3.0144094

Mostrar código

# Añadimos las coordenadas al dataset de entrenamiento para graficar.
train_lda_plot <- train
train_lda_plot$LD1 <- lda_train_pred$x[, 1]
train_lda_plot$LD2 <- lda_train_pred$x[, 2]

head(train_lda_plot[, c("Wine", "LD1", "LD2")])

Mostrar código

# Representación de las observaciones en el espacio discriminante.
ggplot(train_lda_plot, aes(x = LD1, y = LD2, color = Wine)) +
  geom_point(size = 2.5, alpha = 0.85) +
  theme_minimal() +
  labs(
    title = "LDA: proyección de los vinos en el espacio discriminante",
    subtitle = "LD1 y LD2 son combinaciones lineales de las variables químicas",
    x = "Primera función discriminante: LD1",
    y = "Segunda función discriminante: LD2"
  )

Mostrar código

# Gráfico base de MASS. Es útil para visualizar la separación por clases.
plot(wine_lda, main = "Representación de los vinos en el espacio factorial LDA")

Mostrar código

# ==============================================================================
# Proyección de train sobre los ejes discriminantes en Python
# ==============================================================================

X_train_lda_py = lda_py.transform(X_train_py)

lda_train_df_py = pd.DataFrame(X_train_lda_py, columns=["LD1", "LD2"])
lda_train_df_py["Wine"] = y_train_py.reset_index(drop=True)

print(lda_train_df_py.head())

        LD1       LD2 Wine
0 -2.946595  1.296510    1
1 -3.277589  2.139971    1
2  5.170203  2.211306    3
3 -0.514019 -1.710295    2
4 -5.215803  1.942322    1

6.3 Histogramas discriminantes

Los histogramas discriminantes muestran cómo se distribuyen las clases sobre cada eje LD. Son útiles para interpretar si una función discriminante separa bien las categorías.

R
Python

Mostrar código

# ==============================================================================
# Histogramas de las funciones discriminantes
# ==============================================================================

# LD1 suele capturar la mayor parte de la separación.
ldahist(lda_train_pred$x[, 1], g = train$Wine, ymax = 1)

Mostrar código

# LD2 captura separación adicional entre clases.
ldahist(lda_train_pred$x[, 2], g = train$Wine, ymax = 1)

Mostrar código

# En Python construiríamos histogramas por clase usando matplotlib.
# Se deja el código listo para ejecutarse en un entorno Jupyter/Quarto.

import matplotlib.pyplot as plt

for clase in lda_train_df_py["Wine"].cat.categories:
    subset = lda_train_df_py[lda_train_df_py["Wine"] == clase]
    plt.hist(subset["LD1"], alpha=0.5, label=f"Clase {clase}")

plt.title("Histograma de LD1 por clase")
plt.xlabel("LD1")
plt.ylabel("Frecuencia")
plt.legend()
plt.show()

Mostrar código

for clase in lda_train_df_py["Wine"].cat.categories:
    subset = lda_train_df_py[lda_train_df_py["Wine"] == clase]
    plt.hist(subset["LD2"], alpha=0.5, label=f"Clase {clase}")

plt.title("Histograma de LD2 por clase")
plt.xlabel("LD2")
plt.ylabel("Frecuencia")
plt.legend()
plt.show()

6.4 Clasificación con `predict()`

Aunque el script original muestra una clasificación manual basada en reglas sobre LD1 y LD2, en la práctica habitual se utiliza predict(). Esta función devuelve:

La clase predicha.
Las probabilidades posteriores de pertenecer a cada clase.
Las coordenadas de cada observación en el espacio discriminante.

R
Python

Mostrar código

# ==============================================================================
# Predicción sobre test
# ==============================================================================

lda_test_pred <- predict(wine_lda, newdata = test[, -1])

# Clase predicha para las primeras observaciones.
head(lda_test_pred$class)

[1] 1 1 1 1 1 1
Levels: 1 2 3

Mostrar código

# Probabilidades posteriores para las primeras observaciones.
head(lda_test_pred$posterior)

           1            2            3
4  1.0000000 1.431268e-12 2.387577e-17
11 1.0000000 1.794398e-10 2.315713e-19
16 1.0000000 9.353097e-09 7.602262e-13
18 0.9999925 7.473744e-06 5.224970e-13
22 0.9972565 2.743287e-03 2.131227e-07
34 0.9999999 9.907016e-08 2.100790e-15

Mostrar código

# Coordenadas LD1 y LD2 para test.
head(lda_test_pred$x)

         LD1      LD2
4  -4.522271 3.708218
11 -5.064289 2.008654
16 -3.194724 2.783339
18 -3.201553 1.127780
22 -1.576776 1.190986
34 -3.905780 1.535292

Mostrar código

# ==============================================================================
# Matriz de confusión y accuracy para LDA
# ==============================================================================

MC_lda <- table(
  Real = test$Wine,
  Predicho = lda_test_pred$class
)

MC_lda

    Predicho
Real  1  2  3
   1 11  0  0
   2  0 14  0
   3  0  1  8

Mostrar código

accuracy_lda <- sum(diag(MC_lda)) / sum(MC_lda)
accuracy_lda

[1] 0.9705882

Mostrar código

# Versión con caret.
confusionMatrix(
  data = lda_test_pred$class,
  reference = test$Wine
)

Confusion Matrix and Statistics

          Reference
Prediction  1  2  3
         1 11  0  0
         2  0 14  1
         3  0  0  8

Overall Statistics
                                          
               Accuracy : 0.9706          
                 95% CI : (0.8467, 0.9993)
    No Information Rate : 0.4118          
    P-Value [Acc > NIR] : 3.92e-12        
                                          
                  Kappa : 0.9548          
                                          
 Mcnemar's Test P-Value : NA              

Statistics by Class:

                     Class: 1 Class: 2 Class: 3
Sensitivity            1.0000   1.0000   0.8889
Specificity            1.0000   0.9500   1.0000
Pos Pred Value         1.0000   0.9333   1.0000
Neg Pred Value         1.0000   1.0000   0.9615
Prevalence             0.3235   0.4118   0.2647
Detection Rate         0.3235   0.4118   0.2353
Detection Prevalence   0.3235   0.4412   0.2353
Balanced Accuracy      1.0000   0.9750   0.9444

Mostrar código

from sklearn.metrics import confusion_matrix, accuracy_score, classification_report

# ==============================================================================
# Predicción sobre test con LDA en Python
# ==============================================================================

y_pred_lda_py = lda_py.predict(X_test_py)
y_prob_lda_py = lda_py.predict_proba(X_test_py)

print("Matriz de confusión LDA:")

Matriz de confusión LDA:

Mostrar código

print(confusion_matrix(y_test_py, y_pred_lda_py))

[[12  0  0]
 [ 0 13  1]
 [ 0  0 10]]

Mostrar código

print("Accuracy LDA:", accuracy_score(y_test_py, y_pred_lda_py))

Accuracy LDA: 0.9722222222222222

Mostrar código

print(classification_report(y_test_py, y_pred_lda_py))

              precision    recall  f1-score   support

           1       1.00      1.00      1.00        12
           2       1.00      0.93      0.96        14
           3       0.91      1.00      0.95        10

    accuracy                           0.97        36
   macro avg       0.97      0.98      0.97        36
weighted avg       0.97      0.97      0.97        36

6.5 Clasificación manual basada en reglas LDA

El fichero LDA.R incluía una regla manual construida a partir de LD1 y LD2. Esta idea es didácticamente interesante porque ayuda a entender que LDA transforma los datos a un espacio donde las clases pueden separarse mejor.

En este ejemplo se mantiene esa lógica, pero debe entenderse como una regla aproximada dependiente de la base de datos y de la partición train/test.

Mostrar código

# ==============================================================================
# Regla manual basada en LD1 y LD2
# ==============================================================================

# Añadimos LD1 y LD2 al conjunto test.
test_regla <- test
test_regla$LD1 <- lda_test_pred$x[, 1]
test_regla$LD2 <- lda_test_pred$x[, 2]

# Regla inspirada en el script original:
# - Si LD2 < 0, asignamos clase 2.
# - Si LD2 >= 0 y LD1 < 0, asignamos clase 1.
# - En caso contrario, asignamos clase 3.
#
# IMPORTANTE: esta regla no es universal. Funciona solo si la geometría del
# espacio discriminante tiene una estructura compatible con estos cortes.

test_regla$Regla <- NA

for (i in 1:nrow(test_regla)) {
  if (test_regla[i, "LD2"] < 0) {
    test_regla[i, "Regla"] <- "2"
  } else {
    if (test_regla[i, "LD1"] < 0) {
      test_regla[i, "Regla"] <- "1"
    } else {
      test_regla[i, "Regla"] <- "3"
    }
  }
}

test_regla$Regla <- factor(test_regla$Regla, levels = levels(test$Wine))

MC_regla <- table(
  Real = test_regla$Wine,
  Predicho = test_regla$Regla
)

MC_regla

    Predicho
Real  1  2  3
   1 10  1  0
   2  0 14  0
   3  0  2  7

Mostrar código

accuracy_regla <- sum(diag(MC_regla)) / sum(MC_regla)
accuracy_regla

[1] 0.9117647

6.6 Visualización de fronteras de decisión LDA

Para representar fronteras de decisión en dos dimensiones, seleccionamos unas pocas variables. No podemos dibujar una frontera en 13 dimensiones directamente, por eso esta visualización es una simplificación didáctica.

Mostrar código

# ==============================================================================
# Fronteras de decisión LDA con algunas variables
# ==============================================================================

# partimat() es útil para docencia porque muestra regiones de clasificación.
# Solo se ejecuta si el paquete klaR está instalado.

if (requireNamespace("klaR", quietly = TRUE)) {
  partimat(
    Wine ~ Alcohol + Acl + Mg + Proline,
    data = wine,
    method = "lda",
    prec = 200,
    image.colors = c("darkgoldenrod1", "snow2", "skyblue2"),
    col.mean = "firebrick"
  )
} else {
  message("El paquete 'klaR' no está instalado. Se omite partimat().")
}

7 Comprobación de supuestos

LDA y QDA son métodos basados en supuestos probabilísticos. En la práctica, no siempre se cumplen perfectamente, pero conviene conocerlos.

7.1 Supuestos principales de LDA

Las variables predictoras son cuantitativas.
Las clases son conocidas en entrenamiento.
Dentro de cada clase, las variables siguen aproximadamente una normal multivariante.
Las matrices de covarianzas de los grupos son aproximadamente iguales.

7.2 Test de igualdad de matrices de covarianzas

El test de Box’s M contrasta si las matrices de covarianzas de los grupos pueden considerarse iguales. Es especialmente relevante para decidir entre LDA y QDA.

Mostrar código

# ==============================================================================
# Test de Box's M
# ==============================================================================

# H0: las matrices de covarianzas de los grupos son iguales.
# H1: al menos una matriz de covarianzas difiere.
#
# Si el p-valor es pequeño, hay evidencia contra la igualdad de covarianzas.
# Esto puede sugerir que QDA es más adecuado que LDA, siempre que haya datos
# suficientes para estimar las matrices por grupo.

if (requireNamespace("biotools", quietly = TRUE)) {
  biotools::boxM(train[, -1], train$Wine)
} else {
  message("El paquete 'biotools' no está instalado. Se omite boxM().")
}

8 QDA

8.1 Ajuste del modelo QDA

QDA se ajusta también con el paquete MASS, usando la función qda(). La sintaxis es muy parecida a lda(), pero el modelo estimado es distinto porque permite una matriz de covarianzas diferente para cada clase.

R
Python

Mostrar código

# ==============================================================================
# Ajuste del modelo QDA
# ==============================================================================

wine_qda <- qda(Wine ~ ., data = train)

wine_qda

Call:
qda(Wine ~ ., data = train)

Prior probabilities of groups:
        1         2         3 
0.3333333 0.3958333 0.2708333 

Group means:
   Alcohol Malic.acid      Ash      Acl        Mg  Phenols Flavanoids
1 13.75771   2.022917 2.457083 16.82708 105.39583 2.833542  2.9783333
2 12.28614   1.914211 2.247368 20.44035  92.35088 2.272105  2.1182456
3 13.17974   3.434615 2.450769 21.44872  98.25641 1.704615  0.7966667
  Nonflavanoid.phenols  Proanth Color.int       Hue       OD   Proline
1            0.2872917 1.904583  5.455000 1.0562500 3.171458 1115.3542
2            0.3721053 1.604035  3.138070 1.0497544 2.775439  512.5965
3            0.4441026 1.149487  7.352308 0.6815385 1.701795  618.8462

Mostrar código

# ==============================================================================
# Elementos principales del objeto QDA
# ==============================================================================

# Probabilidades a priori de cada clase.
wine_qda$prior

        1         2         3 
0.3333333 0.3958333 0.2708333

Mostrar código

# Medias de grupo.
wine_qda$means

   Alcohol Malic.acid      Ash      Acl        Mg  Phenols Flavanoids
1 13.75771   2.022917 2.457083 16.82708 105.39583 2.833542  2.9783333
2 12.28614   1.914211 2.247368 20.44035  92.35088 2.272105  2.1182456
3 13.17974   3.434615 2.450769 21.44872  98.25641 1.704615  0.7966667
  Nonflavanoid.phenols  Proanth Color.int       Hue       OD   Proline
1            0.2872917 1.904583  5.455000 1.0562500 3.171458 1115.3542
2            0.3721053 1.604035  3.138070 1.0497544 2.775439  512.5965
3            0.4441026 1.149487  7.352308 0.6815385 1.701795  618.8462

Mostrar código

from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis

# ==============================================================================
# Ajuste del modelo QDA en Python
# ==============================================================================

qda_py = QuadraticDiscriminantAnalysis()
qda_py.fit(X_train_py, y_train_py)

QuadraticDiscriminantAnalysis()

Mostrar código

print("Clases:", qda_py.classes_)

Clases: [1 2 3]

Mostrar código

print("Probabilidades a priori:", qda_py.priors_)

Probabilidades a priori: [0.33098592 0.40140845 0.26760563]

Mostrar código

# ==============================================================================
# Predicción con QDA sobre test
# ==============================================================================

qda_test_pred <- predict(wine_qda, newdata = test[, -1])

# Clase predicha para las primeras observaciones.
head(qda_test_pred$class)

[1] 1 1 1 1 1 1
Levels: 1 2 3

Mostrar código

# Probabilidades posteriores.
head(qda_test_pred$posterior)

           1            2             3
4  1.0000000 5.379554e-19 1.548453e-156
11 1.0000000 3.306167e-14 2.740470e-146
16 1.0000000 9.241121e-12 5.002207e-102
18 0.9999999 7.067087e-08 5.075207e-117
22 0.9998904 1.096035e-04  6.823951e-65
34 1.0000000 1.569565e-09  1.024028e-88

Mostrar código

# ==============================================================================
# Matriz de confusión y accuracy para QDA
# ==============================================================================

MC_qda <- table(
  Real = test$Wine,
  Predicho = qda_test_pred$class
)

MC_qda

    Predicho
Real  1  2  3
   1 11  0  0
   2  0 14  0
   3  0  0  9

Mostrar código

accuracy_qda <- sum(diag(MC_qda)) / sum(MC_qda)
accuracy_qda

[1] 1

Mostrar código

confusionMatrix(
  data = qda_test_pred$class,
  reference = test$Wine
)

Confusion Matrix and Statistics

          Reference
Prediction  1  2  3
         1 11  0  0
         2  0 14  0
         3  0  0  9

Overall Statistics
                                     
               Accuracy : 1          
                 95% CI : (0.8972, 1)
    No Information Rate : 0.4118     
    P-Value [Acc > NIR] : 7.908e-14  
                                     
                  Kappa : 1          
                                     
 Mcnemar's Test P-Value : NA         

Statistics by Class:

                     Class: 1 Class: 2 Class: 3
Sensitivity            1.0000   1.0000   1.0000
Specificity            1.0000   1.0000   1.0000
Pos Pred Value         1.0000   1.0000   1.0000
Neg Pred Value         1.0000   1.0000   1.0000
Prevalence             0.3235   0.4118   0.2647
Detection Rate         0.3235   0.4118   0.2647
Detection Prevalence   0.3235   0.4118   0.2647
Balanced Accuracy      1.0000   1.0000   1.0000

Mostrar código

# ==============================================================================
# Predicción sobre test con QDA en Python
# ==============================================================================

y_pred_qda_py = qda_py.predict(X_test_py)
y_prob_qda_py = qda_py.predict_proba(X_test_py)

print("Matriz de confusión QDA:")

Matriz de confusión QDA:

Mostrar código

print(confusion_matrix(y_test_py, y_pred_qda_py))

[[12  0  0]
 [ 0 14  0]
 [ 0  1  9]]

Mostrar código

print("Accuracy QDA:", accuracy_score(y_test_py, y_pred_qda_py))

Accuracy QDA: 0.9722222222222222

Mostrar código

print(classification_report(y_test_py, y_pred_qda_py))

              precision    recall  f1-score   support

           1       1.00      1.00      1.00        12
           2       0.93      1.00      0.97        14
           3       1.00      0.90      0.95        10

    accuracy                           0.97        36
   macro avg       0.98      0.97      0.97        36
weighted avg       0.97      0.97      0.97        36

8.3 Visualización de fronteras de decisión QDA

Igual que antes, visualizamos fronteras usando un subconjunto de variables. En QDA las fronteras pueden curvarse porque cada clase tiene su propia matriz de covarianzas.

Mostrar código

# ==============================================================================
# Fronteras de decisión QDA con algunas variables
# ==============================================================================

if (requireNamespace("klaR", quietly = TRUE)) {
  partimat(
    Wine ~ Alcohol + Acl + Mg + Proline,
    data = wine,
    method = "qda",
    prec = 200,
    image.colors = c("darkgoldenrod1", "snow2", "skyblue2"),
    col.mean = "firebrick"
  )
} else {
  message("El paquete 'klaR' no está instalado. Se omite partimat().")
}

9 Comparación final entre LDA y QDA

R
Python

Mostrar código

# ==============================================================================
# Comparación de rendimiento entre LDA y QDA
# ==============================================================================

resultados <- data.frame(
  Modelo = c("LDA", "QDA", "Regla manual sobre LD1/LD2"),
  Accuracy = c(accuracy_lda, accuracy_qda, accuracy_regla),
  Error = c(1 - accuracy_lda, 1 - accuracy_qda, 1 - accuracy_regla)
)

resultados

Mostrar código

ggplot(resultados, aes(x = Modelo, y = Accuracy, fill = Modelo)) +
  geom_col(alpha = 0.85) +
  ylim(0, 1) +
  theme_minimal() +
  labs(
    title = "Comparación de accuracy entre modelos",
    x = "Modelo",
    y = "Accuracy en test"
  ) +
  theme(legend.position = "none")

Mostrar código

resultados_py = pd.DataFrame({
    "Modelo": ["LDA", "QDA"],
    "Accuracy": [
        accuracy_score(y_test_py, y_pred_lda_py),
        accuracy_score(y_test_py, y_pred_qda_py)
    ]
})

resultados_py["Error"] = 1 - resultados_py["Accuracy"]
print(resultados_py)

  Modelo  Accuracy     Error
0    LDA  0.972222  0.027778
1    QDA  0.972222  0.027778

10 Interpretación docente de los resultados

En este problema, la base de datos Wine suele presentar una separación bastante clara entre las tres clases. Por eso LDA puede funcionar muy bien: las funciones discriminantes capturan direcciones lineales que separan eficazmente los grupos.

QDA puede obtener un rendimiento igual o incluso superior si las clases tienen dispersión distinta. Sin embargo, al estimar una matriz de covarianzas para cada clase, QDA necesita más información y puede ser menos estable en bases de datos pequeñas o con muchas variables.

Desde una perspectiva docente, conviene remarcar:

LDA es más interpretable porque ofrece funciones discriminantes lineales.
QDA es más flexible porque permite fronteras no lineales.
La elección entre LDA y QDA no debe hacerse solo por accuracy, sino también por estabilidad, interpretabilidad y validez de los supuestos.
En problemas con muchas variables y pocas observaciones, QDA puede fallar o sobreajustar.
Si las covarianzas entre clases son parecidas, LDA suele ser preferible.

11 Ejercicios propuestos

11.1 Ejercicio 1

Repite el ajuste de LDA usando solo las variables:

Alcohol
Flavanoids
Color.int
Proline

Compara el accuracy con el modelo que utiliza todas las variables.

11.2 Ejercicio 2

Ajusta QDA con el mismo subconjunto de variables del ejercicio anterior. Compara:

Matriz de confusión de LDA.
Matriz de confusión de QDA.
Accuracy de LDA.
Accuracy de QDA.

11.3 Ejercicio 3

Cambia la semilla de la partición train/test y observa si los resultados se mantienen estables.

11.4 Ejercicio 4

Interpreta los coeficientes de wine_lda$scaling. ¿Qué variables parecen contribuir más a la separación de las clases?

11.5 Ejercicio 5

Estandariza las variables antes de aplicar LDA y compara los coeficientes. Reflexiona sobre por qué la escala afecta a la interpretación de los pesos discriminantes.

12 Código de cierre reproducible

Este bloque resume el flujo mínimo necesario para entrenar y evaluar LDA y QDA en R.

Mostrar código

library(MASS)
library(caret)

url <- "https://gist.githubusercontent.com/tijptjik/9408623/raw/b237fa5848349a14a14e5d4107dc7897c21951f5/wine.csv"
wine <- read.csv(url, sep = ",", stringsAsFactors = TRUE)
wine$Wine <- as.factor(wine$Wine)

set.seed(123)
trainIndex <- createDataPartition(wine$Wine, p = 0.8, list = FALSE)
train <- wine[trainIndex, ]
test <- wine[-trainIndex, ]

modelo_lda <- lda(Wine ~ ., data = train)
pred_lda <- predict(modelo_lda, newdata = test[, -1])
table(test$Wine, pred_lda$class)

modelo_qda <- qda(Wine ~ ., data = train)
pred_qda <- predict(modelo_qda, newdata = test[, -1])
table(test$Wine, pred_qda$class)

	solver	'svd'
	shrinkage	None
	priors	None
	n_components	None
	store_covariance	False
	tol	0.0001
	covariance_estimator	None

	priors	None
	reg_param	0.0
	store_covariance	False
	tol	0.0001

1 Definición del problema

1.1 Contexto

1.2 Objetivo

2 Fundamentos teóricos

2.1 Idea general del análisis discriminante

2.2 LDA: Linear Discriminant Analysis

2.3 QDA: Quadratic Discriminant Analysis

2.4 Diferencias entre LDA y QDA

3 Carga y preparación de datos

4 Análisis exploratorio inicial

5 Partición en entrenamiento y test

6 LDA

6.1 Ajuste del modelo LDA

6.2 Proyección en el espacio discriminante

6.3 Histogramas discriminantes

6.4 Clasificación con predict()

6.5 Clasificación manual basada en reglas LDA

6.6 Visualización de fronteras de decisión LDA

7 Comprobación de supuestos

7.1 Supuestos principales de LDA

7.2 Test de igualdad de matrices de covarianzas

8 QDA

8.1 Ajuste del modelo QDA

8.2 Predicción y evaluación QDA

8.3 Visualización de fronteras de decisión QDA

9 Comparación final entre LDA y QDA

10 Interpretación docente de los resultados

11 Ejercicios propuestos

11.1 Ejercicio 1

11.2 Ejercicio 2

11.3 Ejercicio 3

11.4 Ejercicio 4

11.5 Ejercicio 5

12 Código de cierre reproducible

6.4 Clasificación con `predict()`